Efficient quantum transport simulation for bulk graphene heterojunctions 
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The quantum transport formalism based on tight-binding models is known to be powerful in dealing with a 
wide range of open physical systems subject to external driving forces but is, at the same time, limited by the 
memory requirement's increasing with the number of atomic sites in the scattering region. Here we demonstrate 
how to achieve an accurate simulation of quantum transport feasible for experimentally sized bulk graphene 
heterojunctions at a strongly reduced computational cost. Without free tuning parameters, we show excellent 
agreement with a recent experiment on Klein backscattering [A. F. Young and P. Kim, Nature Phys. 5, 222 
(2009)]. 
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I. INTRODUCTION 

Electronic transport is one of the important fields among the 
increasing number of fundamental studies 1,2 of graphene, a 
one-atom-thick carbon honeycomb lattice. 3 Due to the gapless 
and chiral nature of its electronic structure, graphene exhibits 
energy dispersions linear in momentum, the transport carriers 
behave like massless Dirac fermions, and the properties based 
on Schrodinger wave mechanics in semiconductor physics 
have to be retreated by Dirac-type physics in graphene. Tun- 
neling across pn and pnp junctions is perhaps the most pop- 
ular example that shows how different the charge carriers be- 
have, compared to semiconductor heterostructures. By solv- 
ing the Dirac equation, perfect transmission at normal inci- 
dence across a potential step 4 as well as a potential barrier 5 
was shown for monolayer graphene. This mimicks the Klein 
paradox in quantum electrodynamics 6 and was later referred 
to as Klein tunneling, 7,8 which attracted both experimental 9 " 16 
and further theoretical 13,17-24 investigations. 

The Dirac theory, an effective approach valid only for low- 
energy excitations, generally serves as a starting point for the- 
oretical studies of transport in graphene and can often provide 
analytical results to capture basic physical insights for cer- 
tain problems with simplified system geometries. For further 
considerations, such as to maintain the lattice information on 
graphene or to account for complicated geometries and more 
realistic factors, one has to resort to more advanced theoret- 
ical models. The tight-binding model (TBM), a commonly 
used semiemperical approach for electronic structure calcula- 
tions in solid state physics, 25 allows for consideration of more 
complete band information on graphene at a low computa- 
tional cost. The combination of the TBM with nonequilib- 
rium Green's function approaches forms the modern quantum 
transport formalism, 26 which is able to deal with a wide range 
of conductors composed of a scattering region and external 
leads with or without bias. The description of the graphene 
scattering region of interest, however, requires a TBM Hamil- 
tonian matrix, 
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H gm (V,t,t') = J^V„clc„-t c\c m -t' c t c m, (1) 
n=l (m,n) {{ m , n )) 



sites /V and therefore imposes a computational limit when ad- 
dressing realistic experimental system sizes. This is partly the 
reason why many quantum transport studies address graphene 
"nanoribbons" rather than large-area graphene. The notation 
in Eq. (1) is described as follows: t (f') is the nearest (next 
nearest) neighbor hopping parameter, V„ is the local potential 
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whose matrix size depends on the involved number of atomic 



FIG. 1. (Color online) (a) Schematic of double-gated graphene. 
(b) Carrier density profile n(x) (top) and its corresponding local 
Fermi level Ef(x) (middle). The extracted potential profile V(x) 
(bottom) is given by the difference between the global Fermi level 
Ep and Ep(x); see text, (c) Reproduced densities n(x) provided 
in the Supplementary Material for Ref. 15 with V(, g = 50 V and 
V tg = —8.9, —7.9, • ■ ■ ,0.1, 1.1 V (curves from bottom to top), and the 
extracted corresponding V(x) (curves from top to bottom). 
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energy at site n, c* n (c„) creates (annihilates) a charge carrier 
at the rath site, and the summation Yt(m,n) (L((m«») mns over 
all m and n site indices that are nearest (next nearest) to each 
other within the scattering region. 

Typical sizes of graphene flakes for experimental transport 
investigations amount to a few microns by a few microns, 
but even a 1 /im x 1 fim graphene flake contains roughly 10 7 
atoms, leading to a spinless single-orbital TBM Hamiltonian 
matrix of more than 10 14 elements that requires an exceed- 
ing memory and hence an unreasonable computation burden. 
TBM-based quantum transport for bulk materials therefore 
requires further improvements to overcome the issue of the 
limited scattering region size. In this paper, we demonstrate 
how an accurate TBM-based transport calculation for bulk 
graphene heterojunctions can be performed without free pa- 
rameters, circumventing the problem of large system scales. 

To achieve such a TBM bulk transport simulation, two 
crucial concepts are required, namely, extraction of a re- 
alistic potential profile and description of a bulk graphene 
scattering region, which are described in Sec. II, where a 
brief summary of the quantum transport formalism is also in- 
cluded (Sec. II C). In Sec. Ill, we revisit and simulate the re- 
cent Klein backscattering experiment 15 for transport through 
double-gated graphene [as depicted in Fig. 1(a)] to compare 
with and to demonstrate our approach. Section IV summa- 
rizes the present work. 



profile. Note that the above model makes use of the linear den- 
sity of states that is normally valid in the experimental range 
of the carrier density, although the energy dispersion based on 
the TBM covers the full range. The energy range beyond the 
Dirac model with a nonlinear density of states can, in prin- 
ciple, be treated within the TBM similarly to the process in- 
troduced above, but this would be relevant only far from the 
energy range of interest. 

A realistic carrier density profile depends on the experimen- 
tal geometry and dielectric material of the gate fabrication. In 
the experiment in Ref. 15, n(x) was obtained from an electro- 
static simulation and empirically described by 

, x / 12.8V tg \ 
"(*)= T-Ti 2.5 +% C bg , (4) 

where 12.8 accounts for the effectiveness of the top-gate rel- 
ative to the back-gate, Cb g ~ 7.23 x 10 10 cnT 2 / V is the clas- 
sical (electron number) capacitance of a 290nm-thick SiC>2 
substrate, and the effective half width of the top-gate is w = 
46 nm. 15 Figure 1(c) shows various carrier density profiles de- 
scribed by Eq. (4), subject to Vb g = 50 V and various V tg , and 
the extracted potential profiles, Eqs. (2) and (3). 



n. THEORETICAL FORMULATION 

A. Extraction of a realistic potential profile 

A theoretical study of transport in graphene, whether based 
on Dirac theory or the TBM formalism, requires the potential 
V(x) as an input, which actually means the local energy offset 
of the Dirac point and is often regarded directly as the elec- 
tric potential. In fact, the application of a gate voltage V g does 
not directly raise the Dirac cone by — eV g (— e being the elec- 
tron charge) but enhances or depletes the carrier density, hence 
raising or lowering the local Fermi level. For double-gated 
graphene [Fig. 1(a)], the combination of a top-gate voltage V tg 
and a back-gate voltage Vbg results in a carrier density pro- 
file n(x) such as that shown in the upper panel in Fig. 1(b). Its 
energy dependence, n(E) = sgn(E)E 2 /[n(h\'f ) 2 ], is obtained 
by integrating the density of states over energy. Defining the 
local Fermi level as 

E F (x) = sgn[n(x)]hv F y/ n\n(x)\, (2) 

one obtains the spatially varying height of the filled states, as 
depicted in the middle panel in Fig. 1(b). In a transport calcu- 
lation, the global Fermi level E F is a fixed quantity. Hence to 
account for the profiles of Ep(x) and n(x), one shifts the local 
band offset by applying a local potential 

V(x)=E F -E F (x), (3) 

as depicted in the lower panel in Fig. 1(b). This completes 
the extraction of the potential profile from the carrier density 



B. Bulk graphene scattering region 

In band theory, the electronic structure of a crystal lattice 
can be solved by applying the Bloch theorem, which allows 
us to reduce the problem with infinitely repeated unit cells to 
only one due to translation invariance along each space di- 
mension. For transport calculations, however, the scattering 
region of interest is composed of a certain finite-size area and 
is generally not translationally invariant. For a large flake of 
double-gated graphene, such as that sketched in Fig. 1(a), the 
transverse dimension (along y) is typically a few microns in 
width so that the edges are of minor importance, and we can 
then assume translational invariance in the y direction. 

Consider bulk graphene oriented with zigzag carbon chains 
along the x direction. Up to nearest neighbor hopping, the 
minimal unit cell can be chosen as one hexagon row, i.e., a 
graphene nanoribbon with zigzag chain number N- = 2 with 
transverse periodicity W = 3a, a w 1.42 A being the bond 
length. The wave function at the bottom site (x,ys\(p) of the 
unit cell is related to that at the top site (x,y F \(p) through the 
Bloch theorem as 28 (x,yr + a\<p) = e ^ {x,ya\(p), implying 
\x,yr){x,yT +a\ = e lk ? w \x,y F ) (x,ys\, where k y is the Bloch 
momentum defined within k y W G [— 71, Tc]. This means that 
a kinetic hopping across the upper boundary of the unit cell 
\x,yr) (x,yr +<z| can be equivalently expressed as a periodic 
hopping \x,yr)(x,yB\ modulated by the phase e' k ? w arising 
from the Bloch theorem. Similarly, one can obtain for the 
lower boundary \x,y B )(x,y B — a\ = e~' k y w \x,yB){x,yr\- Incor- 
porating these periodic hopping terms, the TBM Hamiltonian 
for a bulk graphene scattering region can therefore be written 
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as 

H ba}k (V,f,ky)=H gB r{V,t,0) 

+ (-te ik ' w l t 4 m c Bm +U.o), (5) 
V m / 

where c T (cB m ) creates (annihilates) a charge carrier at the 
top (bottom) edge site of the rath hexagon along x, and 
Hg m (V,t,Q), given in Eq. (1), describes an N z — 2 graphene 
nanoribbon. Note that the above description for a bulk scat- 
tering region is restricted neither to nearest neighbor hopping 
(f' = 0) nor to the material graphene. For the present bulk 
transport simulation, however, next-nearest-neighbor hopping 
does not play an important role and we adopt t = 3 eV and 
t' = throughout Sec. III. 



the spin degeneracy is neglected here, while the valley degen- 
eracy is inherently incorporated in //bulk- 

The retarded Green's function Gr of the scattering region 
at energy E in Eq. (7) is obtained from 

G«(£;fc v ) = £ _ [//buik(y ^ )+EL + ^ r (8) 

where //buik(V,f ;k y ) has been given in Eq. (5) and £/, (Lr) 
is the self-energy due to the left (right) lead composed of 
a semi-infinite repetition of unit cells. Adopting a Schur- 
decomposition-based algorithm for the singular hopping ma- 
trix type, 28 the periodic hoppings as used in //bulk can also be 
included in £/, and Lr, enabling us to study pure bulk-to-bulk 
transmission. The spectral matrix functions T/, with / = L,R, 
in Eq. (7) are given by T/ = ;'(£/ — Ej). 



C. Quantum transport formalism 

The quantum transport simulation in the present work is 
restricted to the linear response regime at zero temperature. 
Thus the Landauer conductance 

8(4) = ^ t T(E° F ;ky)dk y (6) 

LKp J-kf 

is the main object and is obtained by integrating the transmis- 
sion function 

T(E;ky)=Tr(T R G R T L Gl), (7) 

which is equivalent to the Fisher- Lee relation. 29 The Fermi 
wave vector in Eq. (6) is approximated from the low-energy 
linear dispersion by kp =Ep/(hvp) = E F /(3ta/2). Note that 



III. KLEIN BACKSCATTERING EXPERIMENT 
REVISITED 

A. Gate-voltage dependence 

Now we revisit the experiment in Ref. 15 by considering 
the extracted realistic potential V(x) and applying the bulk 
TBM transport formalism introduced above. As shown in Fig. 
1(c), the potential profile saturates at roughly ±200 nm, so 
we consider a scattering region described by H], u \k(V(x),t \k y ) 
with length L x = 400 nm. The transport is solely supported 
by the states at the global Fermi level, which is set to E F = 
Ep (x = ±200 nm). We first investigate the top-gate voltage 
dependence of the single-mode conductance g. In Fig. 2(a), 
we directly compare the oscillating features of our computed 
g with the experimental data Gyk, 21 choosing the measured 
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GYK(Vtg, V bg = 40V) and G YK (V tg ,V hg = 60V) curves as ex- 
plicit examples. In both cases, the general features of the mea- 
sured oscillating conductance are well captured by our TBM 
calculation. The Dirac point position of the locally-gated re- 
gion corresponds to the conductance dip. To the left of this 
minimum the transport is in the npn regime exhibiting Fabry- 
Perot-type oscillations due to interference of backscattered 
waves between the np and the pn interfaces. To the right of 
the dip, the transport enters the nn 'n regime, where graphene 
becomes much more transparent than for npn, resulting in the 
suppression of the interference and the rise in the conduc- 
tance. This conductance asymmetry 9 ' 14 ' 19 ' 30 is the first indi- 
rect feature of Klein tunneling, which results in the decay of 
the transmission with the incident angle in the np regime 4 and 
hence a lower integrated conductance, although the tunneling 
at normal incidence is perfect. 

The single-mode spin-degenerate conductance g from Eq. 
(6) has a maximum of 2e 2 jh and does not reflect the main ef- 
fect of the back-gate voltage that tunes the global Fermi level 
Ep : the modulation of the number of modes M participating in 
transport. For bulk graphene at low energy, M can be approx- 
imated by 2kp/Ak y with Ak y = 2n/L y , where L y is the width 
of the graphene flake. This gives M(E) = 2L y \E\ j(nhvp ). 
While the calculation considers the bulk transport across the 
locally gated region in graphene, the contact resistance R c be- 
tween the electrodes and graphene is not included. To com- 
pare with the full map of the measured GYK(Vtg; ^bg)> we tem- 
porarily adopt a simple model to account for multiple modes 
and contact resistance: G(£°) = {[M(E^)g(E^)}- 1 +i? c } -1 . 
Assuming an effective width L y = 2 jj.m and a low contact re- 
sistance R c = 0.2ki2, we display the calculated top- and back- 
gate dependencies of G(Ep) in Fig. 2(b), which qualitatively 
agrees with Ref. 15. Note that the quadrants of G(V tg , Vb g ) are 
determined by the dependence of the potential profile on V tg 
and Vbg, and do not significantly change with the temporar- 
ily introduced parameters Ly and R c , on which we place less 
stress in the present work. 



B. Low-field magnetotransport 

Finally, we come to a closer analysis of the low-field mag- 
netotransport. For an incoherent graphene pup junction a per- 
pendicular magnetic field leads to the increase in the magne- 
toresistance due to the bending of the electron trajectories. 4 
When the top-gate is narrow enough, such as that in Ref. 15, 
with a width of about 20nm, a coherent graphene /;>«/;> junction 
can be formed. Shytov et al. 18 proposed a clever way to ex- 
perimentally test the existence of Klein tunneling, making use 
of the sign change of the Klein backscattering phase at a weak 
magnetic field, which in turn results in a half-period shift of 
the Fabry-Perot oscillations. Based on this semiclassical treat- 
ment the low-field magnetotransport experiment in Ref. 15 
was regarded as providing evidence of Klein tunneling. In the 
following we show that our tuning-parameter-free TBM cal- 
culation confirms the semiclassical picture and, again, agrees 
well with the measurement. 

The orbital contribution of the external magnetic field B z 




n 2 (10 12 cm" 2 ) 

FIG. 3. (Color online) (a) Oscillating part of the computed conduc- 
tance G osc (n2,B z ) (see text for definition) as a function of the car- 
rier density of the locally gated region ni = n(x = 0) and the ex- 
ternal magnetic field B z . (b) Comparison of computed 

(josc curves 

[solid (black) curves] at various magnetic field strengths with the ex- 
perimental data from Ref. 15 [blue (gray) dots] [dotted gray (blue) 
curves]. 



perpendicular to the graphene plane is incorporated in the 
TBM calculation through the Peierls substitution, 31 while the 
Zeeman term is neglected since the Zeeman splitting is rather 
small compared to E F ? To maintain the transverse (y) transla- 
tion invariance throughout the whole system while also keep- 
ing the longitudinal (x) translation invariance in the leads, we 
consider the Landau gauge of A = (Q,xB z ,Q) only in the scat- 
tering region. Inside the left and right leads, however, con- 
stant gauge field strengths A y = XlB z and = xrB z must be 
considered, respectively, where xl and xr are the position co- 
ordinates of the left-most and right-most atomic site of the 
scattering region, in order to avoid a discontinuity of the vec- 
tor potential. 

Since the expected phase shift stems from Klein backscat- 
tering between the two interfaces inside the locally gated re- 
gion, the potential tail does not play a crucial role and we 
reduce the scattering region length to L x = 150nm. Follow- 
ing the definition of the oscillating part of the conductance 
given in Ref. 15, we process our data on the single-mode 
conductance g by first computing the odd part of the con- 
ductance, G dd(n2,B z ) =g(n2,B z ) —g(—ti2,B z ), and then sub- 
tracting its mean value to obtain G OS c("2,#z) = G dd{n2,B z ) — 
G dd( n 2,B z ). Here «2 = n(x — 0) [see Eq. (4)] is the car- 
rier density of the locally-gated region. The obtained oscil- 
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lation fringes of G osc {n2,B z ) are shown in Fig. 3(a), which 
is, again, qualitatively consistent with Ref. 15. The sudden 
phase shift, which indicates the presence of perfect transmis- 
sion and corresponds to the half-period shift predicted by Shy- 
tov et aZ., 18 occurs at magnetic field strengths between 0.2T 
and 0.4T and is in excellent agreement with Ref. 15. In Fig. 
3(b), the computed G osc is compared with the experimental 
data Go S ^(«2,B z ) 27 at various magnetic field strengths (both 
with offset for clarity). 



TBM transport through a bulk scattering region is significantly 
reduced. Together with the realistic potential profile extracted 
from the carrier density profile of a graphene pnp junction, 
this method provides a confirmation of the experiment in Ref. 
15 and its semiclassical theoretical interpretation, at a low 
computational cost without using free tuning parameters. The 
quantum transport approach presented here for studying bulk 
properties is suitable not only for graphene but also for other 
materials where the TBM works well. 



IV. SUMMARY 
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